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' A new method for calculating the free energy of an inhomogeneous superconductor is presented. 

This method is based on the quasiclassical limit (or Andreev approximation) of the Bogoliubov- 
» ! . de Gennes (or wave function) formulation of the theory of weakly coupled superconductors. The 

method is applicable to any pure bulk superconductor described by a pair potential with arbitrary 
spatial dependence, in the presence of supercurrents and external magnetic field. We find that 
both the local density of states and the free energy density of an inhomogeneous superconductor 
\Q ' can be expressed in terms of the diagonal resolvent of the corresponding Andreev Hamiltonian, 

resolvent which obeys the so-called Gelfand-Dikii equation. Also, the connection between the well 
known Eilenberger equation for the quasiclassical Green's function and the less known Gelfand-Dikii 
equation for the diagonal resolvent of the Andreev Hamiltonian is established. These results are 
used to construct a general algorithm for calculating the (gauge invariant) gradient expansion of the 
O ' free energy density of an inhomogeneous superconductor at arbitrary temperatures. 
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^3 . I. INTRODUCTION 



The most interesting, and also most difficult, problems in the theory of weak coupling (BCS) superconductivity!!] are 
those in which the pair potential (order parameter) .has both spatial and time dependence. Examples of such problenae 
are the electromagnetic response of supercpiiductorsa, relaxation phenomena and collective modes in superconductors^, 
O | vortex motion in bulk, wjperconductorsuLi, quantum twmcling of vorticesQ, phase slips in quasi-one-dimensional 
superconducting wiresETtHl, fluctuation effects above Tjl3, etc. In principle all these phenomena can be described 
| in the framework-pf the microscopic theoey of BCS superconductivity in one of its formulations based on cither 
. Green's functions^, or functional integralscJ, or the Bogoliubov-de Gennes (BdG) equationst£l, i.e., the wave function 
formulation. Unfortunately, such an approach is impractical due to formidable technical difficulties of solving the 
corresponding microscopic equations. The existence of the two well separated energy scales in the problem, namely the 
Fermi energy Ep and the magnitude of the gap function A makes the problem even more difficult as far as numerical 
calculations are concerned. However, if we are interested only in the low energy (or long wavelength) physics of 
superconductors then the significant difference between these two energy scales allows us to employ the quasiclassical 
limit of the above mentioned microscopic theories. The quasiclassical Green's function methodliZl is probably the most 
efficient method developed so far for solving problems involving inhomogeneous, non-equilibrium superconductors. 
Nevertheless, this method has its own limitations too (besides the fact that it is valid only on sufficiently long length 
and time scales, for example, the complicated and counterintuitive boundary conditions used in this method need to be 
determined from the underlying microscopic theory, which often relies on questionable approximations). Therefore, it 
is highly desirable to develop an effective theory of weak coupling superconductivity which technically is fairly simple 
and at the same time is general enough to allow for a correct description of the above mentioned phenomena. Such 
an effective theory exists only close to the critical temperature T c , where the sufl£f conducting order parameter is 
small, and a time dependent Ginzburg-Landau (^EQGL) theory is well establishedlifrtj. Recently, attempts to develop 
a viable TDGL-thepry valid at all temperaturesE3c3 yielded some promising results but controversy concerning this 
subject persisted E3. 

So far, all derivations of TDGL theories have been done by using Green's functions and functional integrals. 
Although these methods are suitable for describing inhomogeneous superconductors in the presence of impurities, 
supercurrents and electromagnetic fields, they usually resort to uncontrollable approximations during the decoupling 
of the higher order Green's functions. These approximations may lead to unphysical solutions corresponding to states 
which cannot be described by any wave function. In fact, it is known that the Green's function method as is usually 
formulated does not provide a complete dynamical description of the superconducting system and, therefore, it needs 
to be extended by some extra criterion ( diffe rent from a variational principle) in order to eliminate the spurious, 
unphysical solutions from the correct oner 7 !! 28 !. A typical example in this respect is related to the ground state of 
the superfluid He 3 . Starting from the same BCS reduced Hamiltonian, one can use at least two different forms (or, 
equivalently, decoupling schemes) for the second order correlation function which, in general, lead to different ground 
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states and quasiparticle excitation spectrum: (1) Gor'kovj-and GalitskiiS have obtained an isotropic ground state 
and excitation spectrum, whereas (2) Anderson and More£3, whose approach corresponds to a BCS type of second 
order correlation function, have obtained an anisotropic ground state and excitation spectrum. Interestingly, the 
ground state energy corresponding to the isotropic state is lower than the ground state energy of the anisotropic 
state, however, the former does not correspond to any state wave function and therefore must be rejectedEH. Note 
that there exist other examples as well, where the Green's function method can lead to an imphysical ground state 
with energy smaller than the one obtained by solving the corresponding Schrodinger equationE3. 

In view of this fact it is natural to consider the wave function, or BdG, formulation of weak coupling superconduc- 
tivity to develop a TDGL theory. A first step in this respect is to derive an expression for the free energy functional 
for an inhomogeneous superconductor in time independent (stationary) situation. Such a derivation is the subject of 
the present work. In this paper we present a new method for calculating the free energy density of an inhomogeneous 
superconductor by employing the quasiclassical limit of the wave function formulation of the theory of superconduc- 
tivity. The method is applicable to any pure bulk superconductor described by a pair potential with arbitrary spatial 
dependence, in the presence of supercurrents and external magnetic field. We show that neither the eigenvalues nor 
the corresponding eigcnfunctions of the BdG Hamiltonian are needed to calculate the free energy density, which can 
be expressed solely in terms of the diagonal resolvent of the corresponding Andreev Hamiltonian, resolvent which 
obeys the so-called Gelfand-Dikii equationtil. One of the main features of our method is that it provides a rather 
simple and systematic way to derive the (gauge invariant) gradient expansion of the free energy density at arbitrary 
temperatures. 

The BdG method has been applied previously in the literature to study the physical nmperties of inhomogeneous 
superconductors. The first attempt in this respect has been undertaken by the Orsay groupLJ'Ej. They have determined 
by solving the BdG equations in the quasiclassicaL^WKBJ or Andreev) approximation the low energy excitations 
in the core of an isolated vortex. Also, de GennestS has shown that close to T c the BdG equations can be solved 
by employing the Raylight-Schrodinger perturbation theory and as a result one obtains the Ginzburg-Landau (GL) 
equations. Later on, Bardeen et alE3 (BKJT) have made a more systematic and detailed analysis of the structure 
of an isolated vortex core. BKJT assumed a variational form for both the pair potential and the vector.-Botential 
and the variational parameters have been determined by minimizing the corresponding free energy. ClearyEJ applied 
the theory of BKJT in the vicinity of superconducting transition temperature T c and, quite surprisingly, besides the 
expected Ginzburg-Landau terms in the free energy functional he obtained several anomalous terms as well. These 
findings have been received with great interest by the superconductivity community and several authors have tried to 
explain the origin of these anomalous termsE3. As a result of these research efforts it has been found that apart from 
the vortex problem anomalous terms in the free energy density, a l so appear in other problems involving inliomogencous 
superconducting systems, such as the healing-length problcrr£30, the N-S proximity junction problemtia, etc. Soon 
after the original work of Bar-Sagi and Kupert£l, who managed to find analytically a self-consistent solution of 
the BdG equations in the Andreev approximation (i.e., the Andreev equations) by using a model pair potential 
A(z) oc tanh(az), an intense search has been started to discover other-_practically more useful pair potentials which 
are self-consistent solutions of the corresponding Andreev equationsCJCj. In fact the existence of these self-consistent 
pair potentials are related to the supersymmetric property oLthe properly transformed Andreev Hamiltonian (see 
Sec. [V]) where the pair potential has the role of superpotentialc3. In can be shown that whenever the pair of potential 
energies generated by the superpotential are shape invariant the eigenstates of the corresponding supersymmetric 
Hamiltonians can be determined analytically by means of simple harmonic oscillator like operator algebra. Apparently 
this simple but rather important observation has not been recognized in the literature. The problem of anomalous 
terms in the gradient expansion of the free energy density has been reconsidered by Hta and Eilenberger and Jacobscil 
(EJ) by using the exact self-consistent solution of certain inhomogeneous superconducting systems. These authors 
demonstrated that the actual origin of these anomalous terms are related to surface terms and terms originating from 
the possible discontinuities of the pair potential or its derivatives. EJ have also developed a beautiful theory for 
calculating the free energy density of a quasi one-dimensional inhomogeneous superconductor in the clean limit and 
in the absence of supercurrents and magnetic field. Also, in a recent work Waxmaro starting from the Fredholm 
(functional) determinant expression of the free energy of an inhomogeneous superconductor has shown that the later 
can be expressed in terms of the determinant of a finite 4x4 matrix. However, no viable method for calculating this 
determinant has been proposed. 

The paper is organized as follows: We begin with a brief review of the BdG method of superconductivity (Sec. ||). 
Next, we express the free energy of a bulk superconductor in terms of the spectrum of the BdG Hamiltonian and the 
distribution function of the quasiparticles (Sec. III). The quasiclassical (Andreev) approximation and the expression 



of the free energy in this limit are presented in Sec. IV. Next, by using the wave function formulation of the theory 
of superconductivity, we describe two different methods for calculating the free energy density and the local density 
of states of an inhomogeneous superconductor. Both methods are based on expressing the free energy density of the 
superconductor in terms of the diagonal resolvent of the so-called Gelfand-Dikii equation. The first method, which is 
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applicable only in the absence of the magnetic field and for a real pair potential with arbitrary spatial dependence, 
is presented in Sec. 0, while the second method, which is more general and applic able for superconductors in the 
presence of supercurrents and magnetic field, is presented in Sec. VI. Finally, Sec. VII is reserved for conclusions. 
Also, the derivation of both scalar and matrix Gelfand-Dikii equations, which play a key role in our calculations of 
the free energy, are provided in two appendixes. 



II. THE BOGOLIUBOV-DE GENNES EQUATIONS 



The Bogoliubov-de Gennes (BdG), or wave function, formulation of the microscopic theory of weak coupling 
superconductors represents an attractive alternative to the widely used Green's function and functional integral 
methods. The BdG method is conceptually simple, requires only knowledge of elementary quantum mechanics, yet it 
is general and powerful. In what follows we apply this method to evaluate the free energy density of an inhomogencous 
conventional s-wave superconductor. i_ . 

Mainly to establish notations, we begin with a brief review of the basic equations of the BdG methodtLll Consider 
a pure bulk superconductor in the presence of a static magnetic field. The system is described by the effective mean 
field HamiltomanE£l 



H 



err 



d A r 



rl>t (r ) H (r) ^ (r) + A(r) (r) ^{ (r) + A* (r) ^ (r) V T (r) + 



|A(r)| ; 
V 



(2.1) 



where A(r) is the (mean-field) pair potential, V is the Gor'kov contact pairing interaction [i.e., V (r — r') = 
V 5 (r — r')], the field operators ^v(r) and ipl(r) destroy and create, respectively, an electron at position r with 
spin orientation a =|, |, and obey the usual fermionic anticommutation relations 



{W(r),W(r-')} - 0, {^(r),^(r')} = 0, {^(r),^(r')} = <W 6 (r - r 1 ) 
and, finally, the kinetic energy operator, measured with respect to the Fermi energy Ep, is given by 

H (r) = -L(p- e -A) 2 -E F , p = ~ihV , 
I m \ c / 



(2.2) 



(2.3) 



where t he v ector potential A(r) is related to the total magnetic field H(r) through the equation H(r) = V x A(r) 
In Eq. (2.1), and throughout this paper, implicit summation over repeated spin or pseudospin indices is assumed. 
The effective Hamiltonian (2.1) can be diagonalized by using the Bogoliubov transformationst3 



(2.4) 



where i labels a complete set of quantum states in the relevant Hilbert space, the 7 and 7^ are the Bogoliubov 
quasiparticle annihilation and creation operators, respectively, and satisfy the fermionic anticommutation rules 



(2.5) 



The Bogoliubov amplitudes ui and Vi ought to be determined by the conditiop_.that the transformations (2.4) di- 
agonalize 7i c ff; they obey the so-called Bogoliubov-de Gennes (BdG) equations^ which can be written in compact 
form 



H {r) A(r) 
A*(r) -H*(r) 



*i(r) = Ei^i(r) , 



(2.6) 



where = (v,i, Vi) T is a pseudo spinor in particle-hole space. Thus, the pair potential mixes coherently the particle 
and hole states and , as a result, the Bogoliubov quasiparticles have a mixed particle and hole like character. After 
diagonalization (2.1) reads 
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(2.7) 



where the ground state energy is given by 



E g = -2^ E i J d3r h(r)\ 2 + j d 3 r 



Mr) 
V 



(2.8) 



According to the expression ( |2.7| ) our system is equi vale nt to an "ideal gas" of Bogoliubov quasiparticlcs with energies 
Ei, which are th e ei genvalues of the BdG equations (2.6). For an arbitrary pair potential A(r) the eigenvalue problem 
determined by (2.6), subject to suitable boundary conditions, is difficult to solve even numerically. 



III. FREE ENERGY 

By definition, the free energy is given by@ 

T = (H cS )-TS + T H , 

where S is the entropy, and 



T H = J d 3 rF H (r) , F H (r) 



\H{r)-H a \ 

8tt 



(3.1) 



(3.2) 



is the positive magnetic field exclusion energy due to the screening supercurrents induced by the applied field H a . Also, 
we have assumed that the temperature distribution across the system is homogeneous. The notation (...) = Tr {p . . .} 
indicates the average over some statistical ensemble described by the density matrix p. In thermal equilibrium 



exp (-Wrff/T) 



(3.3) 



Tr{exp(-H cff /T)} ' 

One defines the mean occupation number of level "i" corresponding to spin orientation a by 

fia = (lllia) , (3.4) 

and one assumes that there is no magnetic ordering in the system such that both spin orientations are equally likely, 



fi — /it f'< 



n 



(3.5) 



It is known that the entropy of an ideal gas of fermionic (quasieX particles, which is not necessarily in equilibrium, 
can be expressed in terms of the mean occupation numbers fi asc3 



S = -2 £ [/ 4 In /( + (1-/01x1(1-/,; 



(3.6) 



where the factor of two accounts f or the two independent spin orientations. 

Thus, inserting Eqs. ( p/7 3.4,|3~6) into (3.1), the free energy of the system, which is a functional of the pair potential 
A(r), the mean occupation numbers fi, and the vector potential A(r), can be written as 



T [A(r), fi,A(r)] = T = E g + 2 £ Ei fi - 2 T [fi h-fi + (1 - fi) h (1 - /*)] + ?h 



(3.7) 



In thermodynamic equilibrium one requires the free energy to be stationary with respect to A, fi, and A. Hence, 
stationarity with respect to: (i) the pair potential yields the so-called gap equation (i.e., the self consistency condition 
for the pair potential) 



A(r) = F(r)(V T (r)^(r)) = 7^«i(f)<(r) (1 - 2 f t ) 



(3.8) 
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(ii) the mean occupation number of the state i (for either two spin orientations) yields the usual Fermi distribution 
function 



fi = [exp(^/T) + l]- 
and, (iii) the vector potential yields the Maxwell equation 

4tt 



V x (V x A(r)) = —j(r) , 

c 



where the supercurrent density is given by 

e 



j(r) = 2- V \f iU *(r)P Ui (r) + (1 - /,) v t {r)Pv*{r) 



P 



A(r). 



(3.9) 



(3.10) 



(3.11) 



In the absence of the magnetic field, the BdG equations ([2.6|) together with Eqs. (3.8) and (3.E) yield the standard 



BCS result corresponding to a uniform and real pair potential A(r) = A Q ; the eigenstates i are plane wave states |fe) 
and 



E k = 



Uk 



VN 




(3.12) 



2N-V2 



where N a is the normal state density of states (for both spin orientations) at the Fermi level, oj c is a cut-off frequency 
of the order of the Debye frequency, and u> n — 7rT(2n + 1) are fermionic Matsubara frequencies. 

In what follows we will be interested in calculating the free energy T for a spati ally varyi ng pair potential and 
magnetic field which do n ot n ecessarily obey the self consistency equations (3J3) and ( |3. 10 - 3. 11). Fo r the moment, we 
assume that the relation ( |3.9D is valid but, later on, we will relax this condition as well (see Sec. VD| ). So, we consider 
a superconductor in which the quasiparticles are in thermal equilibrium but the pair potential and the magnetic field 
may have an arbitr ary spatial variation. In this cas e th e expression of the free energy can be further simplified. 
Inserting (3J)) into ( |3.7|) , and by taking into account (2J3), after some straightforward algebra one obtains 



T 



-2T 



In 2 cosh ■ 



2T 



|A(r) 
V 



(3.13) 



Apparently, in order to calculate the free energy ( 3.1 3| ) it is necessary to know the spectrum {Ei} of the BdG 
Hamiltonian HedG for a given pair [potential A(r) (and boundary condition). Fortunately, this is not the case as 
several authors have already shownE3ui, albeit in the absence of any magnetic field and by assuming that A(r) 
depends only on a single spatial coordinate. Indeed, by employing the identitylj 



**»(£)- n 



rn— — oo 



7T 2 (2m + If 



the free energy ( 3.13 ) can be recast as 
T = 



2T^> 



n 



+ d 3 r 



|A(r)r 
V 



(3.14) 



(3.15) 



where uj m are fermionic Matsubara frequencies. The formal divergence of the above expression of the free energy can 
be eliminated by subtracting from T (i.e., by measuring T with respect to) the free energy J~n of the corresponding 
normal state. Thus, by denoting 5 J- = T — Tn, we have 



^=-2Tj> n 



E? 



-2T 



i ui m >0 
lnDct 



d A r 



|A(r)r 



TL 2 



BdG 



HI 



V 



d A r 



|A(r)r 
V 



(3.16) 
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where, TL is the BdG Hamiltonian corresponding to the normal state of the system (i.e., A = 0), and {e^} denote the 
spectrum of-flo- 

Waxmanc3 has shown that the infinite Fredholm (functional) determinant in Eq. (3.16), which contains in an 
encapsulated form all the information on the one-particle excitation spectrum of the superconductor, can be expressed, 
at least in the case of a quasi one-dimensional inhomogeneous superconductor and in the absence of the magnetic 
field, in terms of a finite 4x4 matrix M. However, the actual evaluation of this matrix M(x), which transports 
eigenfunctions of 7iBdG from x = to x = L (L is the size of the system in the relevant x direction) is quite 
complicated and analytical results are possible only for layered systems with a piecewise constant pair potential. In 
the work by Eilenberger and Jacobsci the Fredholm determinant is calculated in terms of a function £{x) which obeys 
an integral equation of Volterra type. This method seems to be somewhat simpler than Waxman's and allows for 
analytical results (in the quasiclassical limit) in several nontrivial cases and, furthermore, provides a viable procedure 
to obtain the gradient expansion of the free energy density about its equilibrium value. 

In contrast to both above mentioned methods, which are only applicable when A(r) varies along a given direction, 
in the absence of any external field, and with the Bogoliubov quasiparticles in thermal equilibrium with the super- 
conducting condensate, our method of calculating the free energy of an inhomogeneous superconductor is valid for 
an arbitrary A(r), in the presence of an arbitrary static magnetic field, and it can be generalized for an arbitrary 
distribution function fi of the quasiparticles. Our method is based on the quasiclassical approximation of the BdG 
equations which we describe next. 



IV. QUASICLASSICAL (ANDREEV) APPROXIMATION 

Superconductors are characterized by two different energy scales, namely the Fermi energy E F and the amplitude 
of the pair potential (gap function) A G at zero temperature. The length scales corresponding to these energies are 
the Fermi wavelength X F ~ fcp 1 ~ hv F /E F which gives the mean inter-particle distance in the system, and the 
superconducting coherence length £ D ~ hv F /A Q which determines the spatial extent of the pair correlation. Since 
in conventional superconductors E F 3> A (or Af <C £ ), as long as we are interested only in the low energy (or 
long wavelength) properties of the system it is legitimate to employ the quasiclassical approximation of the theory 
of superconductivity. The BdG equations are valid on atomic scale and therefore the spinor wave functions \&j(r), 
which vary on a length scale set by fcp 1 , contain more information than it is necessary to calculate, for example, the 
free energy and free energy density of an inhomogeneous superconductor. In general, this excess of information is 
eliminated at the end of the calculations by integrating out the irrelevant high energy (of rapidly oscillating) degrees 
of freedom. A more practical approach is, however, to eliminate these irrelevant degrees of freedom right at the 
beginning_of the calculations by replacing the BdG equations by their quasiclassical limit, i.e., the so-called Andreev 
equations^. For this purpose, one writes the spinor wave function as a rapidly oscillating phase factor (which 
changes on atomic length scale) times a slowly varying amplitude (which changes on a length scale set by the coherence 
length), i.eH 

*i(r) w $„(r; u) exp (i k F ur) . (4.1) 

Thus, in the quasiclassical approximation, the quasiparticles are moving along classical trajectories which are straight 
lines determined by the unit vector u and the "impact parameter" rj_ (which gives the distance of the quasiclassical 
trajectory from the origin of the coordinate system); the position vector in (f4.l|) reads 



r = xu + r± , (4-2) 
where the impact parameter vector r± is normal to it. Nevertheless, the motion along the quasiclassical trajectories is 



quantized and the corresponding eigenstates are labeled in (4.1) by the quantum number n. So, in the quasiclassical 
approximation the state i is specified by the quantum numbers (n,ii,r±) and the trace with respect to the original 
states i must be evaluated according to the formula 

Furthermore, we have (for brevity we omit the arguments) 

V 2 *j = (V 2 $„ + 2 i k F mV $ n - k F $„) exp (i k F ur) , (4.4) 



and therefore, by using Eq. (2.3) in zero magnetic field, one obtains 



G 



2m 



(V 2 $„ + 2 i kp ttV $„) exp (i kp ur) « vf w • (p$„) exp (i A;f u r) , 



where we have neglected the term involving the Laplacian of <!>„ (Andreev approximation) because 

V 2 $„ 



fcp UV $ T 



(kp^y 1 « 1 



(4.5) 



(4.6) 



According to the notion of minimal coupling, in finite magnetic field in Eq. (4.5) one needs to replace p with P 
V- (e/c) A. 



Note that condition (4.6) may not hold for a small fraction of the total number of quasiclassical trajectories 
characterized by u oriented almost perpendicular to V$ n . The Andreev approximation also fails in spatial regions 
where the pair potential (and/or its derivatives) has discontinuities, e.g., at interfaces, boundaries etc. These non- 
analyticities in A(r) reflect the fact that in such regions the pair potential changes rapidly on atomic scale. Within 
the quasiclassical approximation this kind of behavior of A(r) can be described by (nonintuitive) effective boundary 
conditions which must be derived starting from the underlying microscopic theory which is valid on atomic scale. It 
seems to be well established by now that if one does not account properly for the possible discontinuities in the pair 
potential (and or its cLe£i*a|tiyes) these can lead to unphysical anomalous terms in the corresponding Ginzburg-Landau 
free energy function aE3' E3cj. 

Finally, inserting 4T ) into the BdG equations (|]^) and by taking into account ( |4.5| ), one arrives at the so-called 
Andreev equation 



H A 1> n (r) 



H(x) A(x;u,r±) 
A*(x;u,r_ L ) -H*(x) 



$ n (x;u,r_L) = E n (u,r ± )<S> n (x;u,r ± ) , 



where 



H = H(x) = v-p u ■ ( p — - A) = —ifiVFd x — VF-u- A(x) 
V c / c 



(4.7) 



(4.8) 



Note that the Andreev equations (4.7, i. S) are effectively one-dimensional; the independent variable is x (the position 



along the quasiclassical trajectory), and the other degrees of freedom (u,r±) enter the equation only as parameters. 
This is a key observation which allows us to treat inhomogeneous superconductors characterized by a pair potential 
with arbitrary spatial dependence. 



In terms of the energy spectrum of the Andreev Hamiltonian H.a the free energy ST can be written as [cf. Eq. ( 3.16 )1 



5T= ~2Tnhv F N / d 2 r 



-2TitTiv f N d 2 r 



"IT ' ^ 

n 

dllfi, \ ^ 
"4V ^ 

ui m >0 



inn 

w,„>0 



In Dot 



d 6 r 



|A(r) 



■Hi 



HI 



d A r 



V 

|A(r)| 2 
V 



(4.9) 



In the above expression of the free energy the Fredholm determinant involves only the quantum states along an 
individual quasiclassical trajectory. 

In what follows we derive a relatively simple formula for calculating the logarithm of the above Fredholm determinant 
and, consequently, the free energy. We begin with the case of an inhomogeneous superconductor in the absence of 
supercurrents and magnetic field, where the Andreev equations can be decoupled and, therefore, the calculations are 
fairly simple. The more complicated case of a superconductor in the presence of the magnetic field and supercurrents 
requires a completely new method for calculating the free energy density. This method is presented in Sec. VI. 



V. SUPERCONDUCTOR IN ZERO MAGNETIC FIELD 



A. Free Energy 



A key step in our derivation of the free energy pJji-an inhomogeneous superconductor in zero magnetic field is the 
observation, due originally to Bar-Sagi and KuperESO, that the square of the Andreev Hamiltonian (4.7,4.8) can be 
diagonalized and, therefore, the corresponding Andreev equations for the spinor wave function <& n decouple into two 
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independent Schrodinger like equations. Indeed, by dropping all the arguments for brevity, and assuming without any 
loss of generality a real pair potential, one can write 

H 2 + A 2 [H, A] \ ( ~Ti 2 v 2 d 2 x + A 2 -ihv F (d x A) 
"a = Hi = ( = | , (5.1) 

-[H,A] H 2 + A 2 J \ ihv F (d x A) -H 2 v 2 d 2 + A 2 



Q A can be brought to diagonal form by employing the unitary transformation 

V2 1 




= U- 1 = -L ( 1 \ > (5.2) 




Q' A = U- X VL K U = , (5.3) 



where 

H± = H 2 + A 2 ± h v F A' 

•4% 



h 2 v 2 d 2 x + A 2 ±Ti v F (d x A) . (5.4) 



Thus, the spectrum of 7i\ is given by the combined spectra of the two independent one-dimensional Schrodinger like 
operators H±. 

It is worthwhile noticing that the Hamiltonian fl' A is supersymmetric (SUSY) with A(x) playing the role of 
superpotentiallij. In the language of SUSY quantum mechanics, H + and H- correspond to the fermionic and bosonic 
sectors, respectively, and supersymmetry means that the interchange of these two sectors of fl' A leaves the dynamics 
of the system unchanged— r3jhe most useful properties which result from the SUSY of the Hamiltonian tt' A can be 
summarized as followsE3c3& 

1. The Hamiltonians H± can be expressed in terms of the ladder operators 

Q = -hv F d x + A , Q 1 = hv F d x + A, (5.5) 

as 

H+ = QtQ , H- = QQt . (5.6) 

2. The Hamiltonians H± are positive-semidefinite isospectral (up to a zero mode) operators, i.e., 

H± 4>±,n = E 2 X ct>±, n , (5.7) 

where the eigenfunctions </>+,n and </>-.„ are related through 

1 Q^n, (l>+,n = TE^ Q f <t>-, n , \E n \>0. (5.8) 



j^p— — \ En: 

3. The pairing of the eigenstates of H± fails when E n — 0. A zero mode (eigenstate with zero energy) exists 
whenever one of the wave functions 



,o = AAexp (± J dyA{y)^j 



(5.9) 



is normalizable. Since at most one of the above wave functions is normalizable it is clear that one may have 
only one zero mode belonging to the spectrum of either H + or H-. Indeed, assuming, e.g., that </>+, exists, i.e., 
H + <fi+,o — 0, then 
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{4f+,o\H+\ <p+, ) = {<j) + , |<9 + Q| <i>+,o) = \\Q\MW = o 0_, o cx Q^ + . = o, 

and similarly in the opposite case. The necessary condition that one of </>±,o to be normalizable is that A(x) 
has different signs at x — ±00 along the corresponding quasiclassical trajectory. While for conventional s-wave 
superconductors this condition is difficult to be met in zero magnetic fieldcil, in the case of, e.g., unconventional 
d-wave superconductors A(x) can have different signs at the Iwo opposite sides of a quasiclassical trajectory 
which connects two different lobes of the order parameterE3~E3. When the zero mode is absent we say that 
supersymmctry is spontaneously broken and the ground state of Q' A is degenerate (for a given u and r±). When 
the zero mode exists one has a good SUSY and the zero mode is the ground state of £l' A . 

4. Probably the most useful feature of SUSY quantum mechanics is that it allows us to calculate analytically both 
the spectrum and the eigenfunctions of the partner Hamiltonians H± by means of simple algebraic manipulations, 
provided that the partner potentials U±(x;a ) = A 2 {x\a ) ± ft vf A' (x; a Q ) are shape invarian&, i.e., when 
they obey the conditioned 

U+{x-a ) = U-{x;a 1 )+R{a 1 ) , (5.10) 



where a\ is a new set of parameters uniquely determined from the old ones a via the mapping a\ ~ F (a ), 
and the residual term R(a\) is independent of x. A few examples of superpotentials which yield shape in- 
variant potentials U± are: (i) A(x;a a ) cx a t&nh(r) x), (ii) A(x;a ) cx 1 + a a exp(— rjx), (iii) A(x;a a ) cx 
a a / [1 + exp(— i]x)], and (iv) A (x; a D ) cx a a (l + r]x). For all these model pair potentials the eigenstates of the 
Andreev Hamiltonian can be determined analytically by using the machinery of SUSY quantum mechanicsCJ. 
Once the eigenstates of Ha have been determined it i s po ssible to evaluate numerically the value of the pa- 
rameter 77 by imposing the self-consistency condition (3^). Successful calculations along this line have been 
reported by-Bar-Sagi and KuperEJ'LJ for the pair potential (i), by ClintonEj for case (ii), and by Eilenberger 
and Jacobs^ for cases (iii)-(iv). Of course, in principle, it is possible to obtain analytical results for all known 
nontrivial superpotentials (i.e., A(x) in our case) which lead to shape invariant (or factorizable, in the language 
of In feld and HullEa) potentials U± (x) with the possibility of even satisfying the self-consistency (gap) equation 
(3.8). Unfortunately none of these "super" pair potentials correspond to real physical situations and, therefore, 
we will not pursue here this issue in any further details. Nevertheless, it is fair to recognize the potential use- 
fulness of the application of SUSY quantum mechanics in the study of inhomogeneous superconductors within 
the framework of the Andreev approximation, a fact which to our knowledge has not been fully realized so far 
in the literature. 



Before proceeding any further it is useful to introduce new length C and energy £ units via the definitions 

flV~F 



C = 



A 



to 



and 



£ = A n 



(5.11) 



where A D is a suitably chosen constant pair potential, e.g., the equilibrium BCS gap parameter at the considered 
temperature T. In these new units 



H± = H 2 + A 2 ± A' 
= -dl + A 2 ± d T A 



(5.12) 



It is also convenient to measure the free energy in units of N a A 2 . We shall use these units throughout this section. 

In what follows, the fact that H± are supersymmetric will play no special role. The important thing is that H± 
are independent and Schrodinger like. 

Now we introduce the diagonal resolvents R± of the operators H± which will play the central role in our method 
for evaluating the free energy of an inhomogeneous superconductor in zero field. By definition 



Hence 



R±(x;\) = R±(x;\:ii,r±) 



dx R(x\ A) 



(\-H ± y 



E 



(5.13) 



(5.14) 



and a similar relation holds for R corresponding to the reference state described by the Hamiltonian TL . In Eq. (5.14) 
we have used the shorthand notation 
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R = R+ + R- . 

Next, one integrates both sides of Eq. ( |5. 14 ) with respect to the spectral variable A 



dx R{x; A) 



^ln|A-i?„ 2 



const 



(5.15) 



(5.16) 



The integration constant on the RHS of ( 5.1 6| ) can be eliminated by subtracting from this equation the one corre- 
sponding to the reference state. By introducing the notation 



SR 



R — R Q , 



we get 



Finally, by setting A 



dX 



dx SR(x; A) 



In 



\-El 



(5.17) 



(5.18) 



in terms of the diagonal resolvent SR as 



in this last equation, the logarithm of the Fredholm determinant in (4.9) can be written 



InDet 



HI 



El 



—0J m roo 

d\ 



dx SR(x; A) . 



Thus, the free energy (O) becomes 



ST = 2Ttt / d 2 r ± 



~4tT ^ 

ui m >0 



d\ 



dx SR (x; A; it, rj 



3_ 



d A r 



V 



(5.19) 



(5.20) 



where, for clarity, we have listed all the arguments of the diagonal resolvent. 
Now the free energy density 



SF 



d{5T) _ d{5T) 
d 3 r 



dxd 2 r± ' 



as a functional of the inhomogeneous pair potential, can be readily extracted from Eq. ( 5.20 ) 



SF = 5F[A(r) 



dX SR 



u„>0' 



[A(r)r 
V 



(5.21) 



(5.22) 



where (...) = J n ... means averaging over the directions of the quasiclassical trajectories. Note that the only 

difference between the cases, when the pair potential depends only on one coordinate and when it has an arbitrary r 
dependence, is that in the former case the diagonal resolvent does not depend on the impact parameter r± whereas 
in the latter case it does. The above expression of the free energy density does not contain explicitly either the 
eigenvalues or the eigenfunctions of the Andreev Hamiltonian H.a- All the information about the super cond uctor is 
encapsulated in the diagonal resolvent SR which, however, needs to be determined first in order to make (5.22) useful. 

Since R± are the diagonal resolvents of the one-dimensioaal Schrodinger operators H± = — d 2 + U±, with U± = 
A 2 ± A', they obey the so-called Gelfand-Dikii equationEJO 

- 2R±R± +R± +AR 2 ± (U± - X) = 1 . (5.23) 

For c omp leteness a simple derivation of this equation is provided in Appendix |X] (see also Ref. [n]). Equations ( 5.22 ) 
and (5.23) tell us that the free energy density of an inhomogeneous superconductor can be expressed solely in terms 
of the solution of a nonlinear second order ordinary differential equation. Unfortunately the Gelfand-Dikii equation 
cannot be solved analytically for an arbitrary pair potential. However, both the diagonal resolvent and the free energy 
density can be calculated numerically once some appropriate boundary coiid|itions have been specified. Iii-this respect 
our method of calculating SF is similar to the ones considered by WaxmanEJ and Eilenberger and Jacobseil. However, 
while their methods are applicable only to superconductors described by a pair potential which depends on a single 
spatial coordinate and in the absence of supercurrents and magnetic field, our method is valid for pair potentials with 
arbitrary spatial dependence. Another important feature of our approach is that it provides a simple and systematic 
way for obtaining the gradient expansion of SF for an inhomogeneous superconductor with A(r) varying slowly on a 
length scale £ ^> £ Q . 
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B. Gradient Expansion 



For the normal state the pair potential U± = 0, and (5.23) yields 

Ro,± = - ' 



-A 



(5.24) 



For an arbitrary pair potential the general solution of the Gelfand-Dikii equation ( 5.23| ) can be sought as an asymptotic 
series expansion 



oo 

i? ± (x;A) = -J2 R H X ) ( A2 ~ X Y 



(5.25) 



k=0 



where A = A(x) is the pair potential of the inhomogeneous superconductor. For the uniqueness (up to a sign) of this 
expansion see, e.g., Ref. |3l]. Equation (5.25) is the main ingredient in our derivation of the gradient expansion of SF. 



Our strategy is to express first SF in terms of Rr(x), k = 0, 1, . . ., and then to evaluate these expansion coeffic ients 



The latter task can be accomplished in a systematic way by inserting (5.25) into the Gelfand-Dikii equation (5.23) 
and equating the coefficients of the different integer powers of £ = A 2 — A in the resulting equation. Although this 
method can be used to derive a cumbersome analytical expression for the recursion relation obeyed by the coefficients 
R k (x), in practice it is more convenient to carry out the calculations— by employing a computer software which is 
suitable for sophisticated symbolical calculations, such as Mathematical. 

It is easy to see that the first coefficient Rq = 1. Clearly, for the normal state R^ = 1 and the rest of the coefficients 



vanish identically [cf. Eqs. (5.24) and ( 5.25| )], 
k = 1, 2, . . ., one can write 



r: 



= 0, k = 1,2, 



Thus, if one defines 5Rk = RT, + R 



k ' 



and 



5R(x;X) = 



dX5R{x;X) 



1 \ 1 



VA 2 - A v/=A/ 2 



dX 



X>i? fc (z) (A 2 -A) 



(5.26) 



fe=i 



l A l 



^/A 2 -A y/=Xj 2^ 



5 



dX 



oo (A 2 - X) k+ ^ 



2 ,,,„ l _^ ) + .g _i|a_ 



(5.27) 



Inserting (^.27|) into Eq. (|5.22[) yields 



ui m >0 



5F = 4ttT J2 (M " v/^ + A 2 ) + — 



A 2 



E 

Vfe=l L 



2-kT 



E 

o m >0 



K + a 2 ) 



-fc+i 



SR k {x) 
2k-l 



(5.28) 



The first two terms on the RHS of (5.28) give the well known bulk term contribution to the free energy density of the 
superconducting state with respect to the normal state, while the third term gives the actual gradient expansion in 
term of asymptotic power series of the derivatives of the real pair potential A (a;). 

Following the above mentioned strategy for calculating the expansion coefficients R± , we wrote a Mathematica code 
which evaluates analytically, in a systematic fashion, these coefficients. Here we apply our results to calculate the 
gradient expansion of SF up to the fourth order terms, i.e., 



SF w SF + SF 2 + SF 4 



(5.29) 



where 5F Q is given by the first two terms on the RHS of Eq. ( [5.28 ) . Since 5R\ = there is no first order correction 
to SF. In fact one can easily show, based on symmetry arguments, that all odd order contributions to the gradient 
expansion vanishes identically. This does not mean, of course, that all odd order expansion coefficients <5i?2fc+i are 
equal to zero. 
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To calculate 6F2 one needs the following coefficients 



8R 2 = -(A 2 -2AA") 

4 v / 



SR3 = — f20A 2 A' 2 -A" 2 + 2A'A (3) -2AA (4) 
16 V 



(5.30a) 
(5.30b) 



Note that while 6R2 contains only terms of second order in the small parameter £ /£, the coefficient SR3 contains 
both second and fourth order terms as well. None of the higher order coefficients 8Rk contain other second order 
terms in £ a /£. One of the main features of our method is that it can automatically collect all the terms of a given 
order in the various relevant expansion coefficients SRk- Inserting all the second order terms from Eqs. (5.30) into 
Q5.28 ) one obtains 



8F2 



2nT 



E 

w m >0 



A' 2 - 2 A A' 



A 2 A' 



K. + A 2 ) 



2\3/2 



(5.31) 



One can easily see that the second term (proportional to A") on the RHS of Eq. (5.31) can also be expressed in terms 
of A' 2 . Indeed, we have 



A A" 



A 



dA' 



« + A 2 ) 3 / 2 « + A 2 ) 3 / 2 dx 



d 

dx 



A A' 



3 A 2 A' 



v2^/2 



« + A 2 ) 



2\5/2 



(5.32) 



The total derivative on the RHS yields a surface term upon integration with respect to x which, for a bulk supercon- 
ductor with natural boundary conditions, vanishes. In mote, complex superconductivity problems such surface terms 
may lead to anomalous terms in the free energy functionaS. Nevertheless, it is important to notice that it is always 
possible to express the gradient expansion of the free energy density in terms of even powers of the pair potential and 
its derivatives. Another virtue of the computer implementation of our method is that it can automatically perform 
these partial integ rations and return the final result for SF^ in the desired form. 
Thus, Eq. (5.31) can be rewritten as 



SF? = - ( 2ttT 



ui m >0 U 



A' 



A 2 A 
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(5.33) 



(ujl + A^f 2 (u^ + A 2 ) 5 ^ 
Next, we perform the average over u, i.e., the directions of the quasiclassical trajectories; the relevant expression is 



(/(A) A' 2 ) = (/(A) («■ VA) 2 ) = (n inj ) /(A) (ft A) (d j A) 



■ /(A) Aij (diA)(djA) 



/(A) (VA) 2 



(5.34) 



where /(A) is an arbitrary function of the pair potential. This last result clearly depends on dimensionality; in 
rf-dimensions (nirij) = 2 Inserting the above results into (5.32) one obtains 



6F 2 = ^*TN (Hv F ) 2 ^ 







\5/2 



(vAr 



(5.35) 



where we have used the original units. This expression coincides with the well known Werthamer result (Eq. (129) 
in Ref. ^) for a clean superconductor in the absence of supercurrents and magnetic field, obtained by means of 
many-body Green's functions. 

The complexity of calculating the successive terms in the gradient expansion of 5F increases exponentially with 
the order of the term. Nevertheless, by using the computer implementation of our method we were able to compute 
in matter of minutes the fourth order term 6F4. For this purpose one needs to evaluate the expansion coefficients 
5Rk, k — 3, ... ,6 and then filter out all the fourth order terms in the small parameter £ /£ After collecting all these 
terms, we obtain the following expression for the fourth order term in the gradient expansion of 6F 



6F< = ±(2„T 




5 A 2 JL 



(u4 + A 2 ) u/2 (c^ + A^) 



9/2 



A - 



,2^/2 



A'' 



(5.36) 
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This result has been obtained after dropping irrelevant total derivatives in order to express the final result only in 
terms of even powers of the first and second derivatives of A, the only ones which contribute to the fourth order term 
in the gradient expansion. To obtain the final expression for SF4 all we need to do is to average over the directions of 
the quasiclassical trajectories and to restore the original units. 

Apparently, in the derivation of the results prese nted s o far, the particular form of the Fermi-Dirac distribution 
function of the quasiparticles together with formula ( 3.14 ) were crucial. In what follows we show that this is not the 
case and that our method of evaluating the free energy of an inhomogeneous superconductor can be formulated in a 
more general form which is also applicable for a non-equilibrium distribution function fa of the quasiparticles. The 
basic idea is to express the free energy density in terms of the local density of states corresponding to the Andreev 
Hamiltonian (i.e., along an individual quasiclassical trajectory). 



C. Local Density of States 



In this section we present an alternative derivat ion o f the expressions ( 5.28| ) of the free ene rgy density for the case 
of thermal equilibrium without invoking formula (3.14) but rather rewriting the summation ( |4.3| ) over the complete 
set of states i as 



= irfiv F N Jd 2 r ± y^H£\.. = irhv F N J d 2 r ± p(E) dE . . 



(5.37) 



where the density of states (DOS) along the quasiclassical trajectory determined by (u,r±) is given by 

P {E) = p(E;u,r ± ) = - E n (u,r ± )) . (5.38) 



Here {E n } represent the energy spectrum of the Andreev Hamiltonian. Next, let us define the DOS corresponding to 
the SUSY Hamiltonians H± 



p(E) = p + (E)+p.(E) = Y,HE 2 -El) = ±J2 S (E-E n ) = 



P{E) 
2E 



(5.39) 



Thus, from Eqs. (5.39) and (5.14), by employing the formula 



1 



lim 



1 



1 



e— >0+ X ±ie X 



x±iO+ 

the DOS p(E) can be expressed in terms of the diagonal resolvent R = R + + J2_ as 

p(E) = 2Ep(E) = ImV 

7T — ^ 



E 2 -El+iW 



2 E , 
lim Im / dxR\x;E + ie) 

7T e^0+ 



(5.40) 



Now combining equations ( 3.13 ), ( 5.37 ) and ( |5.40 ) lead us to the following expression for the free energy density 



F = -2ittJ^ dE p(E; r) In ^2 cosh 



4T lim Im / E dE In | 2 cosh 



E 

2f J ' r 

E 



(5.41) 



2T 



(R(x;E 2 +ie)) + 



V 



where, by definition, the local DOS is given by 



p(E;r) 



2E 



Im (R (x;E 2 +i0+)) 



(5.42a) 
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or in conventional units [cf. Eq. ( |5.37| )] 

p(E;r) = 2hv F N Elm (R(x;E 2 +i0+)) 



(5.42b) 



The next step is to subtract from (5.41) the free energy density corresp ondin g to the reference normal state and to 
replace SR in the resulting expression by its asymptotic series expansion ( 5.26 ), i.e., 



SR (x;E 2 +ie) 



1 



1 



,VA 2 -E 2 -ie V-E 2 -ie 
Thus, the free energy density SF — F — F^ becomes 



_j oo 



-k-l 



k=l 



SF = SF - 2 T V lim Im 



fc=l 



EdE- 



In (2 cosh^) ,„ , „ 



(A 2 - E 2 - ie) k+ 2 



where, the zeroth order term in the gradient expansion of SF is given by 

1 



5F 



-4 T lim Im / E dE In 2 cosh 



6^0+ 



2TJ V \/A 2 - E 2 



IE 



+ 



V 



(5.43) 



(5.44) 



(5.45) 



By employing contour integration in the complex plane, i t can be shown (s ee Ap pendix |c]) that Eq. (5.45) coincides 
precisely with the hrst two terms on the RHS of Eq. ( 5.28 ). Equation ( 5.44 ) can be further simplified through 
integration by parts 



/>oo 

Jo (1 



EdE 



{l-E 2 -ier 



In 2 cosh ■ 



E 
2T 



Hence 



SF = 5F + J2 



(SR k (x)) 



k=l 



2k -1 e^0+ 



' 2T{2 



lim Im / dE 



1 f 00 

-i TT / dE ~ 

)Jo (1 



tanh 



{\-E 2 -ieY 



k-i 



tanh 



(l-E 2 -ie) 



(5.46) 



(5.47) 



By using complex contour integration, it can be shown that (see Appendix 



Cfc(T) = lim Im / dE 



tanh 



2 T 



(1 - E 2 - ie 



(ci + ir feH 



UJ m >0 



(5.48) 



lim Re / dE 

£->0 + 



tanh 



2 T 



\JE 2 - 1 (l-E 2 



k-l 



k = 2,3, . 



where E = E + ie. Finally, inserting ( |5.48 ) into Eq. Q5.47 ) leads to our previous result ( 5.28) ) and, therefore, to 
Eq. (5.28) which can be also written as 



SF = SF + f^^\ (,i/? A (.r! 



k=2 



2k- 1 



(5.49) 



The coefficients Cfc(T) can be calculated by using their integral representation ( |5.48 ). By employing the identity 
tanh(-E/2T) = 1 — 2 f(E), where f(E) is the Fermi function, one can separate Cfc(T) into a temperature independent 
and a temperature dependent part; the T independent part can be calculated analytically with the result 



Cfe(T) = lim Re / dE 
_ 2 fc - 2 (fc-2)l 



1-2/OE) 



\/E 2 - 1 (l-E 2 



k-l 



(5.50) 



{2k -3)1 



- lim Rc / dE 

£^0 + 



2/(£) 



\JE 2 - 1 (l-E 2 



k-l 
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Furthermore, by repeated partial integration, the second term on the RHS of Eq. (5.50) can be expressed as an 



improper definite integral involving the derivatives of the Fermi function and the familiar BCS DOS 

ME, - ^ . ,5,:, 

For convenience we list below the expressions of the coefficients Cfc (T) for k — 2 and 3 



2 2 
c,(T) = i - J 



l?°^ iE {-%)-\l?° {E)dE {%) 2 • <^ 3) 



where p s {T) is the superfluid density at temperature T. 

As we have already mentioned, this second metho d of ca lculating the free energy of an inhomogeneous superconduc- 



tor by means of the effective local density of states ( 5.42a ) is quite general and in fact it is applicable for an arbitrary 



distribution fa of the Bogoliubov quasiparticles, as we show in the next section. 

D. Non-equilibrium Free Energy Density 

Consider a superconducting state in which the quasiparticles are out of equilibrium with the con den sate. We also 



assume that the superconducting state can be described by the effective mean-field Hamiltonian (2.1), with a pair 



potentia l A(r) and a non-equilibrium quasiparticlc distribution function f(E;r). Then, the local DOS p(E;r) given 



by Eq. (5.42b) is applicable with the same diagonal resolvent R± studied in the previous sections. Thus, one can 



immediately write down the expressions for the energy (W) and entropy (S) densities of the system 

[A(r)l 2 f°° 

W = LJ ^ ZL - / EdE p(E; r) [1-2 /(£;!■)] , (5.54) 



and 



S = dE p(E; r) {f(E; r) In f(E; r) + [1 - f(E; r)} In [1 - f(E; r)]} . (5.55) 
Jo 

The usefulness of these equations depends on the problem at hand. For example, if the system is in local thermal 
equilibrium, such that a local temperature T(r) can be defined and the distribution function of the quasiparticles can 
be expressed, e.g., as f(E; r) — (exp [E/T(r)] + 1) , then it make sense to define a free energy density through the 
usual thermodynamic relation F — W — T S. Furthermore, assuming that the considered superconducting state is 
close to the equilibrium BCS state, it is straightforward to derive a gradient expansion formula for SF along the line 
discussed in the previous sections. 

VI. SUPERCONDUCTOR IN THE PRESENCE OF THE MAGNETIC FIELD AND SUPERCURRENTS 

A. Free Energy 

For A(r) = |A| e t9 complex, with a general spatial dependence of the phase 8, and in the presence of a static 
magnetic field the squared Hamiltonian Ti\ cannot be rotated into a matrix with second-order differential operators 
on the diagonal and off-diagonal terms equal to zero. Consequently we will go back to the expression for the free 



energy (4.£), but now written as 



T = -2Tnhv F N I d 2 r± _ I ^ V InDet (iuj m + H A ) + f d 3 r ^^- + T H . (6.1) 



4tt ^ 'J V 

Here we use the factorization w", + E\ = [iu> m + E n ) (—iu m + E n ), so the sum now goes over both positive and 
negative Matsubara frequencies. 
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The determinant stays unchanged when we make the unitary transformation 



Ha 



-ihv F d x -v F £A u \A\e w 

\A\e-* e iTiv F d x -v F ^A t 



e -i0/2 o 



e 10 / 2 , 
-ihv F d x + v F mv s 



ie/2 q 
e~ lB l 2 



ifivp d x + v F mv s 



(6.2) 



Here 



and 



u ■ v s 



A,. = u ■ A 



u- — V0- — A 

2m \ Tic 



are, respectively, the components of the vector potential and the superfluid velocity in the direction of the quasiclassical 
trajectory. In what follows, it is more convenient to use the Hamiltonian Ha instead of Ha- 

As in the absence of the magnetic field and s uperc urrents, the Fredholm determinant in Eq. ([0]) can be calculated 
from the trace of the matrix resolvent [cf. Eq. ( 5.1 9| ) ] 



G(x, y; r±,u; A) 



We have 



InDct (iuj m + Ha) = ^2 ln tt Um + = ~y~! 



E n — A 



d\ Tr G(A) = - 



dX dxtrG(x,x;r±,u;X) 



(6.3) 



In Eq. (S.3) the notation "Tr" means the trace of the differential operator, and so involves integration over x, while 
the symbol "tr" means a summation over the two spin indices only. In the following, for brevity we will omit the 
arguments rj_ and u. 

As shown in Appendix |bJ the 2x2 matrix G(x, x; A) is obtained from the matrix function g(x; A) = G(x, x; A) 173, 
which satisfies the Eilenberger equation 



ihv F g 



A — v F m v s 
IAI 



IAI 



-A + v F m v Su J ' 



= 0. 



The Eilenberger (or quasiclassical Green's) function g(x; A) also satisfies 

1 



tr.g = 0, g z 
In terms of g(x; A), the free-energy density is given by 



4Ti 2 v 2 F 



ctq 



F = 2Tirhv F N 



■in 



E 



d\tr{g(x; A) <r 3 } 



\Mr)l 
V 



H 



(6.4) 



(6.5) 



(6.6) 
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B. Gradient Expansion 

To obtain the gradient expansion of the free energy density F, we rewrite the Eilenberger equation ( 16.41 ) in the form 

ihv F g'+[V,g] + [A,g] = 0, (6.7) 



where 



is of the zcroth order and 



A = 



A -|A| 
IAI -A 



V = 



-VFrnv Su 
Dp inv s ^ 



is of the first order in gradients of the order parameter. 

The contribution of the zeroth order to g satisfies the equation 

[A, go] = 0, (6.8) 

and so it is at every point the Eilenberger function for a homogeneous superconductor with constant real order 
parameter that would be equal to |A| at that point. This may be calculated explicitly from the resolvent. We find 

9o(x;X) = / , (6.9) 

y y ' 2hv F v/A 2 - |A| 2 Vl A l - A / 



which indeed satisfies Eq. (ji. 



To obtain the higher-order terms, it is convenient to transform to the basis of eigenve ctors of A. That is we find a 
matrix B such that B~ 1 AB = diag. The eigenvalues of A are ±£ with ( = -JX 2 — |A| 2 , and 

B = (\ + J (6-10) 



. |A| A + C. 

A general matrix M = M^a Q + ... + M^a 3 (where a ai a = 1,2, 3, are the Pauli matrices, and Oq is the 2x2 identity 
matrix) transforms into M = M^ao + ... + M^aj, = B~ 1 MB, where 

Af(°) = Af(°) (6.11a) 
M« = M (1) (6.11b) 
(M^\ If A i|A|V^^ (6 . llc) 



[m® ) - I {-i\A\ A J ) 

This transformation is complex-orthogonal rather than unitary, because it was supposed to rotate — i | A| o-i + A 03 
with one component purely imaginary into £03, as it, indeed, does. 
Next we define 

R[x;X) = B- 1 g(x;X)B . (6.12) 

Then 

B- 1 g B = R' +[B- X B\ R] . (6.13) 



From (6.10), we obtain 



D-i T3i A|A|' 

B B - ^ a\ + Bq , 



where B is proportional to the unit matrix and, therefore, contributes nothing to the commutator with R in 
Eq. (6.13). Hence, in this new basis, Eq. (p/fl) reads 
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ihv F R'(x; A) + [U, R(x; A)] + C [cr 3 , R(x; A)] 
where U = B^V B + ihv F (B^ 1 B') {1) is given by 



0, 



f/(3) 



i | A | wf ™ v s 

c 

A v F m v Su 

c 



and R satisfies the conditions 



txR(x; A) ee , ft 



2 _ 



1 



4&V 



^0 



(6.14) 

(6.15a) 
(6.15b) 
(6.15c) 

(6.16) 



It is this matrix function R(x; A) that is the analogue of the scaZajc-resolvent ( 5.13 ); hence we denote it by the same 
letter. The function R can be expanded into the asymptotic serieS 



r — 'Y2 ^ n £ 



n=0 



where 



Ro = B 1 g B 



2Hvf 



<?3 



(6.17) 



(6.18) 



We have somewhat generalized Dikii's workril, because the expansion parameter C" 1 is dependent, and has to 
be, therefore, differentiated too when we substitute ( 5.17] ) into (6.14). The derivative of the n— th order term in the 
expansion is 



(RnC 



R' n C n -nR n ('C 



If we multiply £ by a constant C, bot h terms will be multiplied by C ", so they are both of the order n in £ . 
Equating the n— th order term in ( 6.14 ) to zero, we obtain the recurrence relation 



ih vf I R' r 



c 

n — ft 

c 



[U, R n ] + [a 3 , R n+1 ] = 



(6.19) 



If we write R n = Rn^&o + ■ ■ 
following recurrence relations 



i?i 3 V3, then R^ = since tri? = 0, and the remaining components satisfy the 



R 



(i) _ 

n+l — 



R 



(2) _ 
n+l ~~ 



(ft^C n ) 



= 2 



h vf 



(u^ r^-u^rW) 



(6.20a) 
(6.20b) 
(6.20c) 
„ - n . For 

C, = const, the theory in Ref. |3l] guarantees that the right-hand side of Eq. ( |6.20c ) is always a derivative of a polynomial 
in entries of U and their derivatives. The integration is, therefore, always possible, but it leaves an undetermined 
constant in every R„ \ These constants t ogeth er with the constants in R^ (which are set to zero in our c ase s ince 
tr R = 0) determine the solution of Eq. ( 6.14 ) uniquely. The product of two such solutions again solves ( |l4| ), so 
all the solutions of (6.14) form an infinitely dimensional commutative algebra over the field of complex numbers. 
Equivalently, they form a 2-dimensional algebra over the field of formal series in C^ 1 with constant coefficients. 

For £ spatially dependent, a simple extension of Dikii's theory shows that the right-hand side of (3.20c) can still be 
integrated; now it will contain also powers and derivatives of However, the spatial dependence of £ forces all the 



(3) (3) 

Note that there is no recurrence relation for the coefficients R n , but only for the derivative of Rn C 
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constants on the diagonal of R n to be zero for n > 0, so the solution of (3.14) is completely determined by constants 
on the diagonal of i?o- Moreover, the only solution with Rq = <7o is Co itself, so an arbitrary solution can be written 
as 

R = R 0) a + R { 3) R , 

where R^ and R^ are complex numbers, and R is the unique solution with Rq = 03. Hence, the spatial dependence 
of £ keeps the algebra two-dimensional, but reduces the coefficient field from formal infinite series to complex numbers. 
The algebra is therefore reduced to the two-dimensional Clifford algebra C1(1,C). In our case, R = (i/2h~VF) R- The 
algebra structure then forces R 2 = — (l/47i 2 v F ) ao, in agreement with ( |6.16 ). 

In terms of the expansion coefficients Rf^ the free-energy density has the form 

F = 4T**t* tf £ £ - |A|fl " )( -' A j + t AJ? " )( ' T;A) + ^ + F B . (6.21) 

So, to evaluate the free e nergy density we just need to find t he coefficients i?" from the recurrence relations ( 6.20| ) 
with the initial condition (6.18), substitute them into (6.21), and perform the Matsubara sum and the A— and 
u— integrations. All the A— integrals are of the form 

(V^-IAI 2 

where k is a non-negative integer. 



For k = 0, the integral diverges, and therefore needs special treatment. If we subtract from (6.21) the free-energy 
density of a normal metal Fv then the difference of the corresponding integrals becomes finite 



iX iX 



h, = J dX ^ A2 _ |A|2 - = W - ^» ■ (6.23) 
Using R Q 3 ^ = i/2hvF, we find the zeroth-order contribution to the free energy density 

^ D , v |A|2 

F (r) -F N = 4nTN ^ (cj m - ^ + |A| 2 j + L±- + F H , (6.24) 



ui m >0 



where the spurious divergence of the infinite frequency sum can be eliminated, as usually, by cutting it off at the 
Debye frequency u>d- 

The first-order term -Fi(t*) vanishes because it contains one vector tt to be averaged over the unit sphere which 
gives zero. For k > 1, 



I\,k — 



- K + IAI 2 )-^ . (6.25) 



2k- 

Finally, for the average over the direction it, we use the symmetric integration formula 

, E •«(«s))"-(«(ir, fc _ 1 )-»( 1 r«)) 

("w ' M ) ' ' ' ^ ■ = ^fcTT^ ■ (6 ' 26) 

In this expression many of the terms will be the same. Indeed when vectors vm are all different the numerator on the 
RHS of Eq. ([T2E) contains only {2k)\/k\2 k distinct terms rather than (2k)\. Note that for an odd number of vectors, 
the it-integral vanishes. 

Using Mathematical we obtained the expansion of the free energy density functional up to the eighth order in 
gradients of the order parameter. The terms are getting progressively longer, so we list them below only up to the 
fourth order. To make the formula shorter, we do not perform the it— averaging in the fourth order, just denote it by 
(. . .) around the 21 fourth-order terms. Leaving the explicit directional averaging to the reader has been customary 
in the literature. Also, in the fourth-order terms we write primes instead of gradients. As an example, (v s |A|'Vg) 
means 
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i£ (Mr)(9i|A(r)|)^(r) + T,. i (r)(0 i |A(r)|)ftt; aj -(r) + Us< (r)(^-|A(r)|)9^ si (r)) 



In this notation, the expansion up to the fourth order reads: 



F(r)-F N = AnTN ( w ™ - 



|A|' 



F H +F 2 (r)+F 4 (r) 



where 



F 2 (r) = \vTN m 2 v 2 F £ 



|A| 2 



|A| 2 ) 



3/2 



±nTN h 2 v 2 F J2 



« + |A|" 



^5/2 



(V|A|) 2 



and 



F 4 (r) 



tr\4 (^+|Ap) 7 / 2 



4 m 4 |A| 2 v\ 



1 4 fi 2 m 2 t; 2 |A| 72 

2 C.,2 _i_ I A I2"\5/ 2 



^ 2 n + |A| 2 ) 

4 ft 4 |A|' 4 



|A| 2 ) 7 / 2 96 (w 2 n + | A |2 )9 / 2 



3 

4 fi 2 TO 2 |A|^ s |A|'< 

(^ 2 n + |A| 2 ) 5/2 

4 h 2 m^AI^IAI" 

K + |A| 2 ) 5/2 + 

3 4n 4 |A| 2 |A| //2 
16 (^ + |A| 2 ) 7 / 2 ' 

1 4?j 4 |A| 2 |A|'|A|' 



(^ 2 n + |A| 2 ) 5/2 

35 4ft 2 m 2 |A| 4 t; 2 |Ar 2 

^ 8 (^ n + |A| 2 ) 9 / 2 " 64 ( w 2 n + | A | 2 ) n / 2 

49 4 h 4 |A| 2 |A|' 4 5 4 h 2 m 2 |A| 3 u s |A|' v' a 



25 4 ft 2 m 2 |A| 2 ?; 2 |A| /2 

"8" + | A |2)7/2 

105 4^ 4 |A| 4 |A|' 4 



+ |A| 2 ) 7/2 

:2 _2 | A 13 .,2 



5 4?i 2 m 2 |A| 3 i; 2 |A|' 



1 4 ft 2 to 2 |A| 2 < 2 

+ 4 + |A|2) 5 /2 ~ 4 (cJ 2 i+ | A |2)V2 

3 4n 4 |A||A|' 2 |A|" 77 4n 4 |A| 3 |A| ,2 |A|" 



1G 



K, + |A| 2 ) 



2^/2 



1 4^ 4 |A|'' 



48 

x2 , 



80 {uJ 2 m + \A\ 2 f 2 
1 4a 4 |A|'|A|'' 



1 v F h 2 m 2 |A| 2 t> s t; 

2 (^ + |A| 2 ) 5 / 2 

1 4?i 4 |A||A|"- 



|A| 2 ) 



2^/ 2 



« + |A| 2 ) 



2^/2 



40 



« + |A| 2 ) 



2^5/2 40 



^ 2 „ + |A| 2 ) 



2^/2 



(6.27a) 



(6.27b) 



(6.27c) 



The second order term ( |6.27b| ) is identical with Werthamer's result (Eq. (129) in Ref. |63| ) for a clean superconductor 
in finite magnet ic fiel d, and for v s = this reduces to our previous result |£3g ). In the same limiting case v s = 
the expression ( 6.27c[ ) of the fourth order term gives, up to a total derivative, the same result as (5.36). The fact 
that the two methods we used to calc ulate F 4 are fully independent of one another gives us confidence iu the validity 
of our results. However, the formula ( 6.27c ) apparently disagrees with the result obtained by TewordtE-i Work is in 
progress to locate and understand the difference between these two results and we hope to report our findings in this 
regard in a future publication. 



C. Local Density of States 



As we have seen in Sec. VC, an alternative route for calculating the free energy density of an inhomogeneous 
superconductor is based on the local DOS. The free energy of a bulk superconductor can be written 

T = ln(2cosh|l)+ [ d 3 r^+F H 

Ei>0 ^ J J 

= -2ttT dEp(E) In ^2 cosh + Jd 3 r^- + T H , (6.28) 



where the DOS p{E) in the quasiclassical approximation reads 
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p(E) = - Im 



2hv F N / d 2 r± dx 



47T 



tr G(x,x;E) 



(6.29) 



and E = E + is. Furthermore, Eqs. ( 6.12| ) and ( |6.10 ) allow us to express the trace in ( |6.29 ) in terms of the diagonal 
resolvent R(x; A) 



trG(x, x; E) — tr g(x;E)as 



= tr 



B R(x; E) B^ 1 a 3 



= 2 



tlAlRW^^ + ER^&E) 



(6.30) 



E 2 - |A| 



Hence, by employing the asymptotic expansion ( |6.17 ) the local DOS p(r; E), and the corresponding free energy density 
F(r) can be written, respectively, in the following form 



p(r;E) 



— Tivf N Im 

7T 



dn^^ i\A\RW(x;E) + ERW{x;E) 



47T 



E 

n=0 



E 2 - IAI 2 



n+l 



(6.31) 



and 



F(r) 



\Thv F N Im 



F H (r) . 



2T / 4tt 



n=0 



- IAI 



71+1 



V 



(6.32) 



Similarly to Eq. ( |3.21 ), f irst w e need to determine the rele vant coefficients i?" from the recurrence relations ( 6.20) ) 
with the initial condition ( 3.18 ), then substitute them into ( 6.32 ), and finally perform the E— and u— integrations. 
All the E— integrals are of the form 



DC: 

Je,h — I m J dE In 2 cosh 



E 



iE 



2T 



Re / dE In 2 cosh- 



E 



E 2 - | A | 2 
E 



2k+l 



2T 



E 2 - | A | 2 



2k+l 



(6.33) 



In Appendix ^ we show that the integrals defined by Eqs. ( |6.22 ) and ( |6.33D are related through 



E 7 



A,fe 



—J E ,k, k = 0,1,2,. 

7T 



(6.34) 



By employing this identity, a direct comparison between Eqs. ( 6.21 ) and ( 6.32 ) shows that the two routes to the free 
energy density give the same result. 



VII. CONCLUSIONS 



In this paper we have presented a general method, based on the semiclassical limit of the Bogoliubov-de Gennes 
(or wave function) formulation of the theory of weak coupling superconductivity, for calculating the (gauge invariant) 
free energy density of an inhomogeneous superconductor with a pair potential with arbitrary spatial variation and in 
the presence of supercurrents and magnetic field. We have shown that the free energy density can be expressed in 
terms of the diagonal resolvent of the Andreev Hamiltonian, the semiclassical limit of the BdG Hamiltonian, which 
obey the so-called Gelfand-Dikii equation. Since the solution of the Gelfand-Dikii equation can be easily expressed 
in terms of an asymptotic series, our method is most suitable for obtaining the gradient expansion of the free energy 
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density when the supreconducting order parameter has slow spatial variations on a length scale set by the BCS 
coherence length. To the best of our knowledge, this is the first time when the gradient expansion of the free energy 
of a clean inhomogeneous superconductor, in the general three-dimensional case and in the presence of supercurrents 
and external magnetic field, has been obtained by employing the wave function (BdG) formulation of the theory of 
superconductivity. Our result .for the second order term in the gradient expansion of the free energy density coincides 
with the result of WerthamerE§E3 obtained more than three decades ago by using Green's functions. Howevp. our 
expression of the fourth order term appears to be somewhat different from Tewordt's Green's function resultEa and 
further investigation is needed to establish the origin of this discrepancy. Nevertheless, since in the zero-field case we 
have arrived at the same result for the fourth order term in the gradient expansion of the free energy by using two 
essentially different methods, we are confident in the viability of our approach and results. 

We have also shown that our method for calculating the free energy of an inhomogeneous superconductor is ap- 
plicable for states far from equilibrium characterized by an arbitrary temperature field and quasiparticle distribution 
function. 
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APPENDIX A: THE GELFAND— DIKII EQUATION 



Consider the one-dimensional Schrodinger operator 

H s (x) = -dl + U{x) , 

defined on the interval x € [a, b] (any of a and b may be infinite), and the associated eigenvalue problem 

Hs(x)ip n {x) = E n ip n {x) , 
corresponding to the homogeneous boundary condition 

aip(x) + /3ij}(x) — 0, for x = a, 6 . 
Let us denote by ip a {ipb) the solution of the equation 

Hs(x)i/>(x) = [-dl + U(x)]iP{x) = Ej>(x), 



(Al) 
(A2) 
(A3) 
(A4) 



where E is an arbitrary real number, which obeys the boundary condition (A3) only at x — a (x — b) but not at the 
other end of the interval. Then the Wronskian of ip a and ipb 



W(E) = 1)' a (x)M*)-M*)il>b(x) 

does not depend on x and its only a function of E. Note that W(E) vanishes only for E = E„ 
The Green's function G(x, y; E) associated to Hs is defined through 

H s (x)G(x,y;E) = H s (y)G(x,y;E) = S(x - y) , 



(A5) 



(A6) 



and can be expressed in terms of the Wronskian (A5) as 

1 



G(x,y;E) = (x 



1 



W(E) 



V 

lf s -E 

Q(x - y) ip b (x) rp a (y) + Q(y - x) rp a (x) ipb(y)} , 



(A7) 



where Q(x) is the step function. One can easily check that (A7) obeys Eq. ( A6) and, because by construction satisfies 
the boundary condition (A3), G(x,y;E) is indeed the Green's function associated to Hs- 
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The diagonal resolvent of H$ for a given energy E is defined in terms of the Green's function as 

1 



R{x;E) = { x 



H s -E 



\^ + \[G{x,x + 5 ] E) + G(x + 8,x,E)] = ^'ffiff 



(A8) 



By taking into account Eqs. (AS) and (A4), the first two derivatives of Re = R(x; E) can be written as (for brevity 
we drop the arguments) 



R' E - 



Ip'a fa + Ipa lp' b 

w 



and 



R' E = 2(U-E)R E + 2 



w 



Then, combining Eqs. (A5), ( A1C ) and flASj) we arrive at 

IT' 



(A9) 



(A10) 



(All) 



Finally, inserting (All) into ( AlCQ and after some rearrangements one obtains the desired Gelfand-Dikii equation 

-2R E R E + Re +4R 2 E (U - E) = 1 . (A12) 



APPENDIX B: THE MATRIX GELFAND— DIKII EQUATION 

In this appendix we will study the Andreev Hamiltonian 

Ha = -ta 3 d x + Aaie 1 " 30 . 
and obtain a relation, analogous to the Gelfand-Dikii equation, obeyed by the diagonal part of its resolvent 

1 



G a /3{x,y;E) = 



H A -E 



(Bl) 



(B2) 



Here the indices a and [3 label components in the two-dimensional Nambu space. 

To derive the Gelfand-Dikii analogue we must first write G a p(x,y; E) in a form similar to the expression we used 
earlier for the resolvent G(x,y;E) of the Schrodinger Hamiltonian. Recall that there we had 



G(x,y;E) 



1 



W(E) 



(B3) 



where W(E) = W (ip a , ipj,) = ^Vfc - V't.V'a is the Wronskian of the two solutions ip a , tpb of the homogeneous Schrodinger 
equation Hip at i, = Eip ay i,, and is independent of x. 

We will begin by constructing the resolvent G a p(x,y; E) for the Andreev Hamiltonian on a finite interval [a, b]. 
If required, the limit of an infinite domain can be taken later. To specify the problem completely we must impose 
boundary conditions on the wave functions at a, b in such a way that the Hamiltonian is self-adjoint. This requires 
that the condition for hermiticity, 



Ha^2 



Ha ipi \i>2 



i ifi\ a 3 iii 







(B4) 



be satisfied in manner that treats tpi and ip2 on an equal footing, thus ensuring that the domain of Ha and H A 
coincide. We impose the vanishing condition at each end separately. Thus 



ip* u if>2u - ^*ii>2l =0, x = a,b 



(B5) 
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where u and I refer to the upper and lower components of ip respectively. We disentangle ■01 , V>2 by dividing by tpiu^i 
and find, for example, 



1p2u 
4>2l 



Eq. (B6) requires that all ip in the domain of Ha obey 



for some real angle 8 a . Similarly 



(B6) 



(B7) 



(B8) 



x—b 



These boundary angles 9 a ,b parameterize the family of possible self-adjoint boundary conditions. Physically one may 
think of them as the phases of the order parameter of a superconductor with an infinitely large energy-gap abutting 
the ends of the interval. 

We now notice that, for any two solutions ipi, ip% of the Andreev eigenvalue problem 



H<7 3 d* + Afiie"^) Vi,2 



£tyi, 2 > ( B9 ) 



the quantity 



(BIO) 

is independent of x. To prove this, simply differentiate w and use Eq. flB9|) . We will see that w(jp 1^2) plays the 
same role for the Andreev equation as the Wronskian plays for the Schrodinger equation. For example, using the 
boundary condition we see that ip\a 3 ipi = = ip2 <J 3' , p2- Combining this with the constancy of w(ip\, ^2), we see that 
w must vanish identically if the two solutions ipi, tp2 are proportional to one another. Conversely, if for two solutions 
of Ha^P — Eip we have w(V>i, V^) = at some point (and hence at all points) in the interval, then 

(Bll) 



or 







= 




_ 1pl2 






^u2 


x—a 



(B12) 



The two solutions therefore satisfy the same differential equation with the same initial boundary condition, and so 
must be proportional. We have therefore shown that w^i,^) provides the same test for linear independence as the 
Wronskian, W(ipi,ip 2 )- 

The resolvent G a p(x, y\ E) will have the form 



G af3 (x,y;E) = A%(y)*%(x) , for x < y 
= (y)**(x) , for x > y 



(B13) 



where ^(x), $f„(x) are solutions of the homogeneous equation and satisfy the boundary condition on the left (L) or 
right (R) hand boundary respectively. The jump- condition obtained by integrating 



{-ia 3 d x +Aa 1 e i ^ e -E)G{x,y;E) = a S(x-y) 



across the point x — y is 



i(a 3 ) aa ,[^,(y)A L -^,(y)Af] = 6 a0 



To solve Eq. flB15|) we define W = ^ L a 3 ^ R and use the conditions "f 1 a 3^ 
multiplying Eq. (jB15|) by *&' a L we find 



(B14) 
(B15) 

^^03^ = 0. For example, on 



24 



This collapses to 



In this manner we obtain 



-zWA R = ^ L {y). 



(B16) 



(B17) 



W* 



w 



for x < y 
for x > y . 



(B18) 



Notice that G a p{x, y; E) = G*p a (y, x; E) as befits the resolvent of a self-adjoint operator. 

For the Schrodingcr problem the Gelfand-Dikii equation applies to the diagonal x = y entry in the resolvent. Our 
matrix-valued G a p(x, y; E) is discontinuous at x = y and, as explained in earlier sections, we must define G a p(x, x; E) 
by taking an average of the left and right-hand limits 



G a p(x,x;E) 



(B19) 



W p W* 

It turns out that G Q ^(x, x; E) is not quite the most convenient quantity to work with. Instead we use the matrix 

g a p(x;E) = G afj <{x,x;E) {a 3 ) pl[j . (B20) 

The utility of this modification is related to the coefficient of d x in the Andreev equation being a 3 instead of the 
identity. 

If one takes the square of the matrix g a p, again using ^ L a 3 1 i' L = ty^ R a 3 ^ R = 0, one finds that 



9a8 



1 



W* 



Now 



Similarly 



sap 1 8 a ex 



(B21) 



(B22) 



(B23) 



Provided that £7 is not an eigenvalue we have W ^ 0, so the two column vectors ^ R and are linearly independent 
and together span the two-dimensional vector space at each point x. Consequently these last two equations are telling 
us that 



1 



9a0 ~ --7<> a 3 . 



We also see that the trace of g vanishes 



g aa = W - 4r W*) = . 

2 \W W* J 

We can therefore find three (generally complex) numbers a\, a 2 , a 3 such that 

a\ + a\ + a 2 = 1 



(B24) 



(B25) 



(B26) 



and 



g aa = (aicri + a 2 a 2 + a 3 a 3 ) . 



(B27) 
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After inserting Eq. ( B19| ) into Eq. (B2C) we may now seek the Andreev-equation analogue of the Gelfand-Dikii relation. 
Because the Andreev equation is f irst-order, we have only to differentiate once with respect to x before we are able 
eliminate the d x \V's using Eq. (B9). We immediately find that 



id x g = [g,a 3 (E-Aa ie ^)} 



(B28) 



The equation ( |B28| ) is well-known in the superconductivity literature as a form of the Eilenberger equation. The 

present derivation is much simpler than those usually adduced. In particular the normalization condition Eq. (B24) 

appears automatically and does not have to be introduced by hand. It also has the added advantage of demonstrating 

that the Eilenberger equation should be regarded as the Andreev problem analog of the Gelfand-Dikii equation. 

-l 



Notice that the position dependent matrix g a p(x) is the diagonal part of 
that it commutes with a 3 (H A — E), i.e., 



<7 3 [H A -E 



Equation ( p28|) asserts 



= a 3 (H A -E) ,g 



(B29) 



APPENDIX C: RELATIONSHIP BETWEEN I XK AND J E k 



To prove Eq. ( 6.34 ) we examine the integrals in ( |6.33 ). Again, for k = 0, the integral diverges. By subtracting the 
contribution of the normal metal from the free-energy density, the integral takes the form: 



oo 

f E 
Je,q = Rc / dE In 2 cosh — 



E 



- 1 



E 2 - I A I 2 



(CI) 



This integral still diverges logarithmically. We can formally integrate by parts 

E " 



Je,o — R-c 



{\JE 2 - | A | 2 - £) In 2 cosh 



•IT 



dE 



E 



Re / — tanh — (;/ E 2 - \A\ 2 - E) 



(C2) 



At E = 0, the boundary term vanishes (as e — > 0+, the contribution from the square root is pure imaginary; the second 
term is zero altogether). As E — ► oo, the boundary term goes to — |A| 2 /4T. To do the regularization consistently, we 
have to discard the boundary term. In the remaining integral, we can replace E by E in the argument of tanh, so we 
can write 



Je,o = \rc J ^t^h^E 2 -\A\ 2 -E) 



(C3) 



O+ie 



We note that the integrand becomes complex conjugate upon E — ► E*. Thus, by extending the lower limit of 
integration to — oo + is, we automatically obtain just the real part: 



Je,o — 



1 



dE 



E 



T tanh— (^E 2 -\A\ 2 -E) 



(C4) 



-00 + Z£ 



We can now formally close the contour of integration around the imaginary upper half axis where tanh(i?/2T) has 
simple poles at the (positive) fermionic Matsubara frequencies, and obtain the result: 



Je,o = ^ iy>m — \/u 



|A| 2 ) 



(C5) 



u m >0 



This result together with fl6.23| ) yield Eq. fl6.34| ) for k = 0. 
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For k > 1 everything is straightforward. The integral Je,i does not contribute because the first-order term is zero 
due to symmetry discussed above. For k > 1, Je,r converges. Integrating by parts and extending the contour gives 



Jem = Re 



■ 1 In 2 cosh 7^ 



2k- 1 



E 2 - I A 



2fc-l 



1 /" dE 



tanh 



2k-l 



(C6) 



The integrated out term vanishes. To perform the integral we again close the contour in the upper half plane and 
obtain 



j£ < fe =4(2^1) ? 2m 



— Y 



-™>o (v/-^- | A p 
1 



2fc-l 



2(2fc- 1) 



^ (i) 2fc - 2 V^+l A I 2 



2fc-l 



(C7) 



By comparing this result with (6.25) one can infer that (6.34) holds for any positive integer k 
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